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INTRODUCTION 

A finite-difference viscous-inviscid interaction program has been developed for simulating the 
separated transonic flow about lifting airfoils, including the wake. In contrast to most interac- 
tion programs, this code combines a finite-difference boundary-layer algorithm with the invis- 
cid program. The recently developed finite-difference boundary-layer code efficiently simulates 
attached and reversed compressible boundary-layer and wake flows. New viscous-inviscid in- 
teraction algorithms were also developed to couple the boundary-layer code and the inviscid 
transonic full-potential program (see ref. [1] for details). Transonic cases with shock-induced 
and trailing-edge separation are computed and compared with experimental and Navier-Stokes 
results. 


VISCOUS ALGORITHM 

The compressible boundary-layer equations for the steady, two-dimensional flow of a perfect 
gas are written in £(x),r)(x, y) coordinates* 

PHH & +Vv r )x) + vu v Tj y ] = -PptU (la) 

pc p [tt(r € £, +T ri ri x ) + vT r ,T) y ] = PupsZ x + {xT^y)^ + /*K»7y) 2 U&) 

{pu)i U + (pu)t, Vx + {pv) n Tfy == 0 (lc) 

where the equations are nondimensionalized as outlined in references [2-3]. By using a general X,y 
to £,77 coordinate transformation, a complex similarity transformation is avoided and a solution- 
adaptive grid is easily employed. In this work, the grid height is varied as a function of the 
computed displacement thickness. 

Near and in reversed-flow regions, the boundary-layer equations are solved in the inverse mode 
to avoid singular behavior at the separation point. Here the wall shear stress Tw and the wake 
centerline velocity u wc are used as the inverse forcing functions because they are efficiently imple- 
mented in the following numerical algorithm. 

The boundary-layer equations are solved using an implicit predictor-corrector algorithm. 
Streamwise marching begins by predicting estimates of u, T, p, v, fi , and K . This predic- 
tor step uses first-order £ operators in the x-raomentum and energy equations, but for 
one step produces second-order accurate values. Because the predictor step only needs to 
be first-order accurate, any nonlinear coefficients can be lagged in A second-order ac- 
curate corrector step is then used to calculate improved values of u, T, p, V, fl, and K . 
During the corrector step, the nonlinear coefficients are evaluated using the most recently 
computed flow variables. Overall, second-order accurate solutions are obtained at the cost 
of two scalar bidiagonal and four scalar tridiagonal matrix inversions per streamwise sta- 
tion. 

The algorithm is presented below in operator notation, where E is the shift operator 

(e.g., Ef'u, = tij—i ): 

v e =(1 - £*-')/( AO A f =(£?> - 1)/(A{) 

6 ( =(Ep -Ef')/( 2A?) S ( =(£+i - £-*)/( Af) 

A predietor-step result is denoted by a tilde (e.g., u) and, unless indicated other- 
wise, the indices are J,k. The superscripts 1,11, and III denote flow-dependent opera- 
tions. 
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Predictor Step 

Predict u at the new station using the x-momentum equation 

+ V*8,U) {pv)j-\-9 , {VyGrfV’) — Ptx^^P “I" Vy$ t,{Pj+ 8> Vy (2fl) 


Predict T at the new station using the energy equation 

Pj+s 1 tfxGriT) Vj-^-gi T)y 6 n 7*] = (26) 

Pii£xG%P Vy^i ,{ k j+ 8' Vy^ t]7) Pj-\-e' iVy^r,^) 

Obtain ~p — p/T and integrate the continuity equation for v 

1 + E~ l tx^zipu) + Vx6*{pv) 

V v (pv) = (2c) 

Z 7/w 


The coefficients p and K are then evaluated using the Cebeci turbulence model [4]. 

Corrector Step 

Correct u at the new station using the x-momentum equation 

MZxfi'lu + ri x 6r,u) + pvirjyStju) = — PZx&$p + Vy'^niPVy'^u) (3o) 


Correct T at the new station using the energy equation 

pc P N6r T + rj x 6, T) + v(rj y 6 V T )] = pu£ x 6^p + r} y (kt) v 5 n T ) + p(r) y 6 V uf (36) 


Update p , v y p, and K as before. 


Flow-dependent operators are used near reversed flow. The logic for determining the 
appropriate operators by specifying values of s 1 and s is summarized in figure 1. To 
incorporate the flow-dependent differencing, the operators 6^ , 6 ^ fl , and 6^'" are defined 
as 
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( 4 ) 


When using flow-dependent differencing values at j -f- 1 and j -J- 2 are obtained from a 
previous iterative sweep. Sweeping of the viscous flow is already required by the viscous-inviscid 
iterations. 


The algorithm is modified to operate in the inverse mode by replacing the pressure tern) 
in the x-momentum equation with expressions containing the inverse forcing functions. These 
relations are obtained by applying the x-momentum equation at the wall and wake centerline. 
For example, at a wall ( k = 1) the x-momentum equation yields an equality between the pressure 
term and the viscous term evaluated at the wall. Replacing the pressure term with a difference 
approximation of the viscous term in, for example, the corrector step x-momentum difference 
equation yields 


pu((x < 5 ? u + Tlx dr, u) + pv(T)y 6r , «) = 


02+01 «3 

-2 

ifc 


— T„ 


+ %M/i»?yM) 


( 5 ) 
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A similar adaptation is made in the wake. 

Equation (5) shows that in the inverse mode u% appears in the difference equation at every k 
index. If this term is treated implicitly, a tridiagonal-like matrix with an additional column of 
nonzero coefficients is obtained. This augmented scalar tridiagonal matrix system is efficiently 
solved using an LU decomposition algorithm (see refs. [2-3]). Because the inverse forcing 
functions are implemented implicitly, no iteration is required to obtain the desired T w or 

VISCOUS-INVISCID INTERFACES 

If T w falls below a prescribed value the viscous algorithm converts from the direct to 
the inverse mode from that point on, including the entire wake. When operating in 

the inverse mode Tw or u wc must be updated so the viscous and inviscid pressures con- 
verge. A number of schemes for updating Tw and u W c were studied. The following are 

the fastest and most reliable of those studied and were used to obtain all the presented 

results: 

r n+l _ rj + fc)(p J _ p"+' ) where: ui = 10 (6a) 

Kt' =<C« - P? + ') w. = 2 (6b) 

A transpiration velocity [5] is used to introduce the influence of the viscous region upon the 
inviscid flow. After 6* has been calculated the transpiration velocity is computed and then con- 
verted to the required perturbations of the inviscid contravariant velocities used in the inviscid 
algorithm. This approach allows the use of inviscid grids that are not orthogonal to a body 
surface or wake centerline in interaction codes 

In many cases it is necessary to account for the viscous flow curvature and the pres- 
sure jump that occurs across these curved stream-tubes. Therefore, the method of ac- 
counting for the pressure variation across the shear layers developed by Lock and Firmin 

[6] is incorporated into the present algorithm The reader is referred to reference [3] for 
details. 

RESULTS 

The interaction code was first tested by comparison with Navier-Stokes computations 

[7] of the M 0 0 = 0.720 flow past an 18% biconvex airfoil. The agreement between the 
present results and the free-flight results of Levy is good in terms of both the computed 
pressure distributions (see fig. 2) and separation points. Also shown are experimental data 

[8] and the results of Levy’s calculations which include the tunnel walls. It is apparent 
that the experimental dat\ were affected by the tunnel walls Experimental data and 
computational results are also available for this airfoil at Moq = 0.754. Experimentally, this 
flow can be unsteady unless a trailing-edge splitter plate is installed, which is modeled in 
the computations. As shown in figure 3, there is some shock-position discrepancy between 
the present computations and the experimental data and results of Levy; it is attributed to 
the wind-tunnel wall effects. Otherwise, the present results are in good agreement with the 
Navier-Stokes results. Both computations predict a greater trailing-edge pressure recovery 
than observed experimentally The computed Cj distribution for this case (see fig. 4) 
indicates that shock-induced and trailing-edge separation are predicted. To indicate the 
convergence rate, the Cf n history for these cases is presented in figure 5. The Moo = 0.720 
case has converged by the 75th iteration, whereas because of the strong shock-separated 
boundary-layer interaction, the M 0 0 = 0.754 case requires approximately 160 iterations to 
converge. 

McDevitt [9] also measured M pea k on the 18% biconvex (with trailing-edge splitter plate) 
versus M^. Figure 6 compares the results of the present method to this experimental 
data. The comparison is encouraging, especially considering the tunnel-wall effects on 
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this data. Also, at the higher Mach numbers some of the discrepancy may be due to 
the isentropic inviscid flow assumption. The large difference between the inviscid and 
viscous peak Mach numbers is indicative of the strong viscous-inviscid interaction being 
modeled. 

To verify the treatment of lifting airfoils, the results of computing the Moo = 0.8 
flow about a NACA 0012 at a = 1° are compared with experimental data [10] and 
Navier-Stokes [11] results. Figure 7 shows good agreement between the experimen- 
tal and computed C p distributions. Both computations predict shock-induced separa- 
tion from approximately 56% to 67% of chord and trailing-edge separation from 95% of 
chord. 

As a final case, the results of computing the = 0.73 flow about the RAE 2822 airfoil 
at Ci = 0.803 are presented. Figures 8-10 compare the C p ,Cf\ e , and 6* distributions 
found experimentally [12], those computed by the present method, and those computed by 
Mehta [13] using a Navier-Stokes code. Mehta performed his computations at a = 2 . 79 ° 
(and computed a Ci — 0 . 793 ); the present computations were performed at a = 2 . 81 ° to 
match the measured Ci = 0 . 803 . The present results are in good agreement with both the 
Navier-Stokes and experimental results. The velocity vectors near the trailing edge computed 
by the present method are presented in figure 11. This figure indicates the high resolution 
obtainable. 

These simulations were obtained on 223x31 inviscid and 223x50 viscous grids. For the cases 
presented, the required Cray-XMP CPU time was 7 to 15 sec, and on the average, 0.0006 sec/grid 
point were required to obtam a converged solution In contrast, an optimized version of the 
thin-layer Navier-Stokes code developed by Steger and Pulliam [ll] requires about 0.06 sec/grid 
point to obtain a converged solution 

SUMMARY 

A fast, versatile, direct-inverse finite-difference boundary-layer code has been developed and 
coupled to a transonic full-potential airfoil code via new viscous-inviscid interaction algorithms. 
The developed interaction code has been used to compute non-lifting and lifting separated 
transonic airfoil flows, and the results are in good agreement with experimental data and 
Navier-Stokes computations. Furthermore, the cost of these solutions is less than twice the 
cost of the inviscid solutions and close to two orders-of-magnitude less than Navier-Stokes 
solutions. 


REFERENCES 

1. Holst, T. L., AIAA Journal, Vol. 17, Oct. 1979, pp. 1038-1045. 

2. Van Dalsem, W. R., and Steger, J. L., AIAA Paper No 83-1689. 

3 Van Dalsem W. R., Ph.D. Thesis, Stanford University, Jun. 1984. 

4. Cebeci, T., AIAA Paper No. 70-741. 

5. Lighthill, M. J., Journal of Fluid Mechanics, Vol. 4, 1958, pp. 383-392. 

6. Lock, R. C., and Firmin, M. C. P., RAE Technical Memo. 1900, Apr. 1981. 

7. Levy, L. L., AIAA Journal, Vol. 16, Jun. 1978, pp. 564-572. 

8. McDevitt, J. B., Levy, L. L., and Deiwert, G. S., AIAA Paper No. 75-878. 

9. McDevitt, J. B., NASA TM-78549, Jan. 1979. 

10. Whitfield, J. D., et al., Computers and Fluids, Vol. 8, 1980, pp. 71-99. 

11. Pulliam, T. H., CFD User’s Workshop, Tullahoma, Tenn., Mar. 1984. 

12. Cook, P. H., McDonald, M. A., and Firmin, M. C. P., AGARD AR 138, 1979. 

13. Mehta, U., Second Symposium on Aerodynamic Flows, Jan. 17-20, 1983. 


4 




•PREDICTOR STEP r^T^ 1 

CORRECTOR STEP = rJJ 
Fig. 1 Flow-dependent differencing logic. 


o EXPERIMENT (8] 



Fig. '3 C p distributions for an 18% biconvex at 
A/oc = 0.754, Rtao = 8*10® 
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Fig. 4 Computed C / distribution for an 18% 
biconvex at Moo = 0 754, Rtoo = 8x10® 
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Fig 2 C v distributions for an 18% biconvex at Fig. 5 Minimum C/ history for the Moo — 0.720, 
Moo = 0.720, f?eoo = 11*10®. i*«<x> = 11*10® and Moo = 0.754, Re oo = 8*10® 

solutions 
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o EXPERIMENT [10] 



Fig. 7 C v distributions for a NACA 0012 
airfoil at Moo = 0 800, Re^ = 2 25*10®, 
0 = 1 ° 


o EXPERIMENT [12] 

NAVIER-STOKES (FREE-STREAM) [13] 



Fig. 8 C p distributions for an RAE 2822 airfoil 
at Moo = 0.730, f?«oo = 6 50*10® ,Cj = 0 803 



X 


Fig. 9 C/ | e distributions for an RAE 2822 
airfoil at Moo — 0.730, Rtoo — 6 50*10®, 
Ci = 0.803. 

o EXPERIMENT [12] 

NAVIER-STOKES (FREE-STREAM) [13] 
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Fig. 10 6 * distributions for an RAE 2822 airfoil 

at Moo = 0 T30,R«oo = 6.50*10® ,C( = 0 803 



Fig. 11 Velocity rectors near the trailing 
edge of an RAE 2822 airfoil at Moo = 0.730, 
R»oc — 6.50*10® ,Ci = 0.803. 
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